home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / C / Applications / Python 1.3.3 / Python 133 SRC / Objects / complexobject.c < prev    next >
Text File  |  1996-01-29  |  11KB  |  574 lines
  1. /* Complex object implementation */
  2.  
  3. /* Borrows heavily from floatobject.c */
  4.  
  5. #ifndef WITHOUT_COMPLEX
  6.  
  7. #include "allobjects.h"
  8. #include "modsupport.h"
  9.  
  10. #include <errno.h>
  11. #include "mymath.h"
  12.  
  13. #ifdef i860
  14. /* Cray APP has bogus definition of HUGE_VAL in <math.h> */
  15. #undef HUGE_VAL
  16. #endif
  17.  
  18. #ifdef HUGE_VAL
  19. #define CHECK(x) if (errno != 0) ; \
  20.     else if (-HUGE_VAL <= (x) && (x) <= HUGE_VAL) ; \
  21.     else errno = ERANGE
  22. #else
  23. #define CHECK(x) /* Don't know how to check */
  24. #endif
  25.  
  26. #ifdef HAVE_LIMITS_H
  27. #include <limits.h>
  28. #endif
  29.  
  30. #ifndef LONG_MAX
  31. #define LONG_MAX 0X7FFFFFFFL
  32. #endif
  33.  
  34. #ifndef LONG_MIN
  35. #define LONG_MIN (-LONG_MAX-1)
  36. #endif
  37.  
  38. #ifdef __NeXT__
  39. #ifdef __sparc__
  40. /*
  41.  * This works around a bug in the NS/Sparc 3.3 pre-release
  42.  * limits.h header file.
  43.  * 10-Feb-1995 bwarsaw@cnri.reston.va.us
  44.  */
  45. #undef LONG_MIN
  46. #define LONG_MIN (-LONG_MAX-1)
  47. #endif
  48. #endif
  49.  
  50. #if !defined(__STDC__) && !defined(macintosh)
  51. extern double fmod PROTO((double, double));
  52. extern double pow PROTO((double, double));
  53. #endif
  54.  
  55.  
  56. /* elementary operations on complex numbers */
  57.  
  58. int c_error;
  59. static complex c_1 = {1., 0.};
  60.  
  61. complex c_sum(a,b)
  62.     complex a,b;
  63. {
  64.     complex r;
  65.     r.real = a.real + b.real;
  66.     r.imag = a.imag + b.imag;
  67.     return r;
  68. }
  69.  
  70. complex c_diff(a,b)
  71.     complex a,b;
  72. {
  73.     complex r;
  74.     r.real = a.real - b.real;
  75.     r.imag = a.imag - b.imag;
  76.     return r;
  77. }
  78.  
  79. complex c_neg(a)
  80.     complex a;
  81. {
  82.     complex r;
  83.     r.real = -a.real;
  84.     r.imag = -a.imag;
  85.     return r;
  86. }
  87.  
  88. complex c_prod(a,b)
  89.     complex a,b;
  90. {
  91.     complex r;
  92.     r.real = a.real*b.real - a.imag*b.imag;
  93.     r.imag = a.real*b.imag + a.imag*b.real;
  94.     return r;
  95. }
  96.  
  97. complex c_quot(a,b)
  98.     complex a,b;
  99. {
  100.     complex r;
  101.     double d = b.real*b.real + b.imag*b.imag;
  102.     if (d == 0.)
  103.         c_error = 1;
  104.     r.real = (a.real*b.real + a.imag*b.imag)/d;
  105.     r.imag = (a.imag*b.real - a.real*b.imag)/d;
  106.     return r;
  107. }
  108.  
  109. complex c_pow(a,b)
  110.     complex a,b;
  111. {
  112.     complex r;
  113.     double vabs,len,at,phase;
  114.     if (b.real == 0. && b.imag == 0.) {
  115.         r.real = 1.;
  116.         r.imag = 0.;
  117.     }
  118.     else if (a.real == 0. && a.imag == 0.) {
  119.         if (b.imag != 0. || b.real < 0.)
  120.             c_error = 2;
  121.         r.real = 0.;
  122.         r.imag = 0.;
  123.     }
  124.     else {
  125.         vabs = hypot(a.real,a.imag);
  126.         len = pow(vabs,b.real);
  127.         at = atan2(a.imag, a.real);
  128.         phase = at*b.real;
  129.         if (b.imag != 0.0) {
  130.             len /= exp(at*b.imag);
  131.             phase += b.imag*log(vabs);
  132.         }
  133.         r.real = len*cos(phase);
  134.         r.imag = len*sin(phase);
  135.     }
  136.     return r;
  137. }
  138.  
  139. complex c_powu(x, n)
  140.     complex x;
  141.     long n;
  142. {
  143.     complex r = c_1;
  144.     complex p = x;
  145.     long mask = 1;
  146.     while (mask > 0 && n >= mask) {
  147.         if (n & mask)
  148.             r = c_prod(r,p);
  149.         mask <<= 1;
  150.         p = c_prod(p,p);
  151.     }
  152.     return r;
  153. }
  154.  
  155. complex c_powi(x, n)
  156.     complex x;
  157.     long n;
  158. {
  159.     complex cn;
  160.  
  161.     if (n > 100 || n < -100) {
  162.         cn.real = (double) n;
  163.         cn.imag = 0.;
  164.         return c_pow(x,cn);
  165.     }
  166.     else if (n > 0)
  167.         return c_powu(x,n);
  168.     else
  169.         return c_quot(c_1,c_powu(x,-n));
  170.  
  171. }
  172.  
  173. PyObject *
  174. PyComplex_FromCComplex(complex cval)
  175. {
  176.     register complexobject *op = (complexobject *) malloc(sizeof(complexobject));
  177.     if (op == NULL)
  178.         return err_nomem();
  179.     op->ob_type = &Complextype;
  180.     op->cval = cval;
  181.     NEWREF(op);
  182.     return (object *) op;
  183. }
  184.  
  185. PyObject *
  186. PyComplex_FromDoubles(double real, double imag) {
  187.   complex c;
  188.   c.real = real;
  189.   c.imag = imag;
  190.   return PyComplex_FromCComplex(c);
  191. }
  192.  
  193. double
  194. PyComplex_RealAsDouble(PyObject *op) {
  195.   if (PyComplex_Check(op)) {
  196.     return ((PyComplexObject *)op)->cval.real;
  197.   } else {
  198.     return PyFloat_AsDouble(op);
  199.   }
  200. }
  201.  
  202. double
  203. PyComplex_ImagAsDouble(PyObject *op) {
  204.   if (PyComplex_Check(op)) {
  205.     return ((PyComplexObject *)op)->cval.imag;
  206.   } else {
  207.     return 0.0;
  208.   }
  209. }
  210.  
  211. complex
  212. PyComplex_AsCComplex(PyObject *op) {
  213.     complex cv;
  214.     if (PyComplex_Check(op)) {
  215.         return ((PyComplexObject *)op)->cval;
  216.     } else {
  217.         cv.real = PyFloat_AsDouble(op);
  218.         cv.imag = 0.;
  219.         return cv;
  220.     }   
  221. }
  222.  
  223. static void
  224. complex_dealloc(op)
  225.     object *op;
  226. {
  227.     DEL(op);
  228. }
  229.  
  230.  
  231. void
  232. complex_buf_repr(buf, v)
  233.     char *buf;
  234.     complexobject *v;
  235. {
  236.     if (v->cval.real == 0.)
  237.         sprintf(buf, "%.12gj", v->cval.imag);
  238.     else
  239.         sprintf(buf, "(%.12g%+.12gj)", v->cval.real, v->cval.imag);
  240. }
  241.  
  242. static int
  243. complex_print(v, fp, flags)
  244.     complexobject *v;
  245.     FILE *fp;
  246.     int flags; /* Not used but required by interface */
  247. {
  248.     char buf[100];
  249.     complex_buf_repr(buf, v);
  250.     fputs(buf, fp);
  251.     return 0;
  252. }
  253.  
  254. static object *
  255. complex_repr(v)
  256.     complexobject *v;
  257. {
  258.     char buf[100];
  259.     complex_buf_repr(buf, v);
  260.     return newstringobject(buf);
  261. }
  262.  
  263. static int
  264. complex_compare(v, w)
  265.     complexobject *v, *w;
  266. {
  267. /* Note: "greater" and "smaller" have no meaning for complex numbers,
  268.    but Python requires that they be defined nevertheless. */
  269.     complex i = v->cval;
  270.     complex j = w->cval;
  271.     if (i.real == j.real && i.imag == j.imag)
  272.        return 0;
  273.     else if (i.real != j.real)
  274.        return (i.real < j.real) ? -1 : 1;
  275.     else
  276.        return (i.imag < j.imag) ? -1 : 1;
  277. }
  278.  
  279. static long
  280. complex_hash(v)
  281.     complexobject *v;
  282. {
  283.     double intpart, fractpart;
  284.     int expo;
  285.     long x;
  286.     /* This is designed so that Python numbers with the same
  287.        value hash to the same value, otherwise comparisons
  288.        of mapping keys will turn out weird */
  289.  
  290. #ifdef MPW /* MPW C modf expects pointer to extended as second argument */
  291. {
  292.     extended e;
  293.     fractpart = modf(v->cval.real, &e);
  294.     intpart = e;
  295. }
  296. #else
  297.     fractpart = modf(v->cval.real, &intpart);
  298. #endif
  299.  
  300.     if (fractpart == 0.0) {
  301.         if (intpart > 0x7fffffffL || -intpart > 0x7fffffffL) {
  302.             /* Convert to long int and use its hash... */
  303.             object *w = dnewlongobject(v->cval.real);
  304.             if (w == NULL)
  305.                 return -1;
  306.             x = hashobject(w);
  307.             DECREF(w);
  308.             return x;
  309.         }
  310.         x = (long)intpart;
  311.     }
  312.     else {
  313.         fractpart = frexp(fractpart, &expo);
  314.         fractpart = fractpart*2147483648.0; /* 2**31 */
  315.         x = (long) (intpart + fractpart) ^ expo; /* Rather arbitrary */
  316.     }
  317.     if (x == -1)
  318.         x = -2;
  319.     return x;
  320. }
  321.  
  322. static object *
  323. complex_add(v, w)
  324.     complexobject *v;
  325.     complexobject *w;
  326. {
  327.     return newcomplexobject(c_sum(v->cval,w->cval));
  328. }
  329.  
  330. static object *
  331. complex_sub(v, w)
  332.     complexobject *v;
  333.     complexobject *w;
  334. {
  335.     return newcomplexobject(c_diff(v->cval,w->cval));
  336. }
  337.  
  338. static object *
  339. complex_mul(v, w)
  340.     complexobject *v;
  341.     complexobject *w;
  342. {
  343.     return newcomplexobject(c_prod(v->cval,w->cval));
  344. }
  345.  
  346. static object *
  347. complex_div(v, w)
  348.     complexobject *v;
  349.     complexobject *w;
  350. {
  351.     complex quot;
  352.     c_error = 0;
  353.     quot = c_quot(v->cval,w->cval);
  354.     if (c_error == 1) {
  355.         err_setstr(ZeroDivisionError, "float division");
  356.         return NULL;
  357.     }
  358.     return newcomplexobject(quot);
  359. }
  360.  
  361.  
  362. static object *
  363. complex_pow(v, w, z)
  364.     complexobject *v;
  365.     object *w;
  366.     complexobject *z;
  367. {
  368.     complex p;
  369.     complex exponent;
  370.     long int_exponent;
  371.  
  372.      if ((object *)z!=None) {
  373.         err_setstr(ValueError, "complex modulo");
  374.         return NULL;
  375.     }
  376.  
  377.     c_error = 0;
  378.     exponent = ((complexobject*)w)->cval;
  379.     int_exponent = (long)exponent.real;
  380.     if (exponent.imag == 0. && exponent.real == int_exponent)
  381.         p = c_powi(v->cval,int_exponent);
  382.     else
  383.         p = c_pow(v->cval,exponent);
  384.  
  385.     if (c_error == 2) {
  386.         err_setstr(ValueError, "0.0 to a negative or complex power");
  387.         return NULL;
  388.     }
  389.  
  390.     return newcomplexobject(p);
  391. }
  392.  
  393. static object *
  394. complex_neg(v)
  395.     complexobject *v;
  396. {
  397.     complex neg;
  398.     neg.real = -v->cval.real;
  399.     neg.imag = -v->cval.imag;
  400.     return newcomplexobject(neg);
  401. }
  402.  
  403. static object *
  404. complex_pos(v)
  405.     complexobject *v;
  406. {
  407.     INCREF(v);
  408.     return (object *)v;
  409. }
  410.  
  411. static object *
  412. complex_abs(v)
  413.     complexobject *v;
  414. {
  415.     return newfloatobject(hypot(v->cval.real,v->cval.imag));
  416. }
  417.  
  418. static int
  419. complex_nonzero(v)
  420.     complexobject *v;
  421. {
  422.     return v->cval.real != 0.0 && v->cval.imag != 0.0;
  423. }
  424.  
  425. static int
  426. complex_coerce(pv, pw)
  427.     object **pv;
  428.     object **pw;
  429. {
  430.     complex cval;
  431.     cval.imag = 0.;
  432.     if (is_intobject(*pw)) {
  433.         cval.real = (double)getintvalue(*pw);
  434.         *pw = newcomplexobject(cval);
  435.         INCREF(*pv);
  436.         return 0;
  437.     }
  438.     else if (is_longobject(*pw)) {
  439.         cval.real = dgetlongvalue(*pw);
  440.         *pw = newcomplexobject(cval);
  441.         INCREF(*pv);
  442.         return 0;
  443.     }
  444.     else if (is_floatobject(*pw)) {
  445.         cval.real = getfloatvalue(*pw);
  446.         *pw = newcomplexobject(cval);
  447.         INCREF(*pv);
  448.         return 0;
  449.     }
  450.     return 1; /* Can't do it */
  451. }
  452.  
  453. static object *
  454. complex_int(v)
  455.     object *v;
  456. {
  457.     double x = ((complexobject *)v)->cval.real;
  458.     if (x < 0 ? (x = ceil(x)) < (double)LONG_MIN
  459.               : (x = floor(x)) > (double)LONG_MAX) {
  460.         err_setstr(OverflowError, "float too large to convert");
  461.         return NULL;
  462.     }
  463.     return newintobject((long)x);
  464. }
  465.  
  466. static object *
  467. complex_long(v)
  468.     object *v;
  469. {
  470.     double x = ((complexobject *)v)->cval.real;
  471.     return dnewlongobject(x);
  472. }
  473.  
  474. static object *
  475. complex_float(v)
  476.     object *v;
  477. {
  478.     double x = ((complexobject *)v)->cval.real;
  479.     return newfloatobject(x);
  480. }
  481.  
  482.  
  483. static object *
  484. complex_new(self, args)
  485.     object *self;
  486.     object *args;
  487. {
  488.     complex cval;
  489.  
  490.     cval.imag = 0.;
  491.     if (!PyArg_ParseTuple(args, "d|d", &cval.real, &cval.imag))
  492.         return NULL;
  493.     return newcomplexobject(cval);
  494. }
  495.  
  496. static object *
  497. complex_conjugate(self)
  498.     object *self;
  499. {
  500.     complex c = ((complexobject *)self)->cval;
  501.     c.imag = -c.imag;
  502.     return newcomplexobject(c);
  503. }
  504.  
  505. static PyMethodDef complex_methods[] = {
  506.     {"conjugate",    (PyCFunction)complex_conjugate,    1},
  507.     {NULL,        NULL}        /* sentinel */
  508. };
  509.  
  510.  
  511. static object *
  512. complex_getattr(self, name)
  513.     complexobject *self;
  514.     char *name;
  515. {
  516.     complex cval;
  517.     if (strcmp(name, "real") == 0)
  518.         return (object *)newfloatobject(self->cval.real);
  519.     else if (strcmp(name, "imag") == 0)
  520.         return (object *)newfloatobject(self->cval.imag);
  521.     else if (strcmp(name, "conj") == 0) {
  522.         cval.real = self->cval.real;
  523.         cval.imag = -self->cval.imag;
  524.         return (object *)newcomplexobject(cval);
  525.     }
  526.     return findmethod(complex_methods, (object *)self, name);
  527. }
  528.  
  529. static number_methods complex_as_number = {
  530.     (binaryfunc)complex_add, /*nb_add*/
  531.     (binaryfunc)complex_sub, /*nb_subtract*/
  532.     (binaryfunc)complex_mul, /*nb_multiply*/
  533.     (binaryfunc)complex_div, /*nb_divide*/
  534.     0,        /*nb_remainder*/
  535.     0,        /*nb_divmod*/
  536.     (ternaryfunc)complex_pow, /*nb_power*/
  537.     (unaryfunc)complex_neg, /*nb_negative*/
  538.     (unaryfunc)complex_pos, /*nb_positive*/
  539.     (unaryfunc)complex_abs, /*nb_absolute*/
  540.     (inquiry)complex_nonzero, /*nb_nonzero*/
  541.     0,        /*nb_invert*/
  542.     0,        /*nb_lshift*/
  543.     0,        /*nb_rshift*/
  544.     0,        /*nb_and*/
  545.     0,        /*nb_xor*/
  546.     0,        /*nb_or*/
  547.     (coercion)complex_coerce, /*nb_coerce*/
  548.     (unaryfunc)complex_int, /*nb_int*/
  549.     (unaryfunc)complex_long, /*nb_long*/
  550.     (unaryfunc)complex_float, /*nb_float*/
  551.     0,        /*nb_oct*/
  552.     0,        /*nb_hex*/
  553. };
  554.  
  555. typeobject Complextype = {
  556.     OB_HEAD_INIT(&Typetype)
  557.     0,
  558.     "complex",
  559.     sizeof(complexobject),
  560.     0,
  561.     (destructor)complex_dealloc,    /*tp_dealloc*/
  562.     (printfunc)complex_print,    /*tp_print*/
  563.     (getattrfunc)complex_getattr,    /*tp_getattr*/
  564.     0,                /*tp_setattr*/
  565.     (cmpfunc)complex_compare,    /*tp_compare*/
  566.     (reprfunc)complex_repr,        /*tp_repr*/
  567.     &complex_as_number,            /*tp_as_number*/
  568.     0,                /*tp_as_sequence*/
  569.     0,                /*tp_as_mapping*/
  570.     (hashfunc)complex_hash,     /*tp_hash*/
  571. };
  572.  
  573. #endif
  574.